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Abstract 


Porcine deltacoronavirus (PDCoV) is a newly emerged swine enteropathogenic coronavirus 
affecting pigs of all ages and causing diarrhea problems. Research findings indicate that 
PDCoV has evolved strategies to escape innate immune response in host cells, but mecha- 
nism of PDCoV in innate immune modulation is not well understood. In this study, we report 
our findings on identifying the alterations of host cell innate immune response affected by 
PDCoV infection and exploring the gene expression profiles of PK-15 cells at 0, 24, and 36h 
PDCoV post infection by RNA sequencing. A total of 3,762 and 560 differentially expressed 
genes (DEGs) were screened by comparison of uninfected PK-15 cells and infected PK-15 
cells at 24 h post infection (hpi) (INF_24h versus NC), and also comparison of infected PK- 
15 cells between 24 and 36 hpi (INF_36h versus INF_24h), which included 156 and 23 por- 
cine innate immune-related genes in the DEGs of INF_24h versus NC and INF_36h versus 
INF_24h, respectively. Gene Ontology function classification and Kyoto Encyclopedia of 
Genes and Genomes signaling pathway enrichment analysis were performed based on the 
DEGs that exhibited the same expression tendencies with most of the innate immune-asso- 
ciated genes among these PK-15 cell samples described above. The enrichment results 
indicated that extensive gene functions and signaling pathways including innate immune- 
associated functions and pathways were affected by PDCoV infection. Particularly, 4 of 5 
innate immune signaling pathways, which were primarily affected by PDCoV, played impor- 
tant roles in /-/FNs antiviral function in innate immune response. Additionally, 16 of the host 
cell endogenous miRNAs were predicted as potential contributors to the modulation of 
innate immune response affected by PDCovV. Our research findings indicated that the innate 
immune-associated genes and signaling pathways in PK-15 cells could be modified by the 
infection of PDCoV, which provides a fundamental foundation for further studies to better 
understand the mechanism of PDCoV infections, so as to effectively control and prevent 
PDCoV-induced swine diarrheal disease outbreaks. 
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Introduction 


Porcine deltacoronavirus (PDCoV) is a newly emerged swine enteropathogenic coronavirus 
that can cause acute diarrhea, vomiting, dehydration, and mortality in pigs [1, 2]. PDCoV was 
first discovered in the rectal swab of pigs in Hong Kong and was initially reported in 2012 [2, 
3]. In 2014, it appeared in the state of Ohio in the United States, resulted in outbreaks of diar- 
rheal disease of piglets and sows accompanied by 30~40% mortality rates in piglets, and later 
rapidly spreading to several other states adjacent to Ohio [4, 5]. At present, PDCoV has been 
identified in multiple countries including Canada, South Korea, Thailand, and China [6-9]. 
Although a variety of PDCoV detection approaches have been established and effectively used 
for diagnostic tests, such as one-step quantitative real-time PCR (qRT-PCR) assays and 
enzyme linked immunosorbent assay (ELISA) [10, 11]. However, there are still no effective 
treatments and control measures for swine PDCoV infections due to its mechanism of 
PDCoV infections is mostly remained unknown. 

Innate immune response provides an essential early line of host cell defense against viral 
infection. When a virus invades a host cell, viral pathogen-associated molecular patterns 
(PAMPs) are recognized by host cell pattern recognition receptors (PRRs), primarily consist- 
ing of toll-like receptors (TLRs), retinoic acid-inducible gene (RIG)-I-like receptors (RLRs), 
nucleotide-binding and oligomerization domain (NOD)-like receptors (NLRs), and cyto- 
plasmic DNA sensors. Subsequently, PRRs stimulate signaling that triggers an innate immune 
response inducing the production of type I interferons (I-IFNs), mainly consisting of several 
IFN-a subtypes and IFN-f. I-IFNs are the major antiviral cytokines of innate immune response 
that mediate related responses resulting in the production of IFN-stimulated genes (ISGs) in 
order to control virus replication [12, 13, 14]. 

To successfully proliferate in host cells, multiple types of coronaviruses have evolved differ- 
ent strategies in order to evade cell innate immune response either by inhibiting the produc- 
tion of I-IFN or by disrupting I-IFN-involved innate immune signaling pathways. Severe acute 
respiratory syndrome coronavirus (SARS-CoV) can induce the formation of acute respiratory 
distress in human beings and has a mortality rate close to 10%. SARS-CoV has gained wide- 
spread attention since it first emerged as a huge threat to humans in 2002. In order to explore 
the pathogenic mechanism of SARS-CoV, various studies have focused on cellular mecha- 
nisms, including innate immune response evasion. Kopecky-Bromberg SA et al. found that 
the nucleocapsid (N) protein which encoded by SARS-CoV againsts IFN-B production by 
inhibiting IRF-3, the key protein necessary for the expression of IFN-£ [15]. Furtherly, Hu Y 
et al. found that N protein of SARS-CoV interfers with RIG-I activation and ubiquitination 
which mediated by tripartite motif protein 25 (TRIM25), resulting in the inhibition of IFN-B 
production which achieve higher SARS-CoV replication in host cells [16]. In addition to 
SARS-CoV, the N protein of other coronaviruses, such as mouse hepatitis virus (MHV), has 
also been identified as an innate immunity antagonist. RIG-I and melanoma differentiation 
gene 5 (MDAS) are both PRRs that can recognize PAMPs and then activate signaling to pro- 
mote I-IFN production. MHV N protein inhibits the synthesis of IFN-B by interacting with the 
protein activator of protein kinase R (PACT), which can bind to RIG-I and MDAS, promoting 
the production of IFN-8 [17]. As an avian coronavirus, infectious bronchitis virus (IBV) 
escapes innate immune response by cleaving the mitochondrial antiviral signaling protein 
(MAVS), which disrupts downstream of RIG-I and MDAS, antagonizing IFN-f synthesis [18]. 

In terms of PDCoV, there are also several strategies have been evolved to escape innate 
immune response. PDCoV acts against the innate immune response by impeding the activa- 
tion of transcription factors IRF3 and NF-«B, both of which are involved in the RIG-I signaling 
pathway disrupting IFN-8 production [19]. Nonstructural protein 5 (NSP5), which is encoded 
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by PDCoV, antagonizes I-IFN signaling through modulation of the Jak/STAT signaling path- 
way by cleaving STAT2, resulting in the disruption of ISGs production [20]. The NF-«B essen- 
tial modulator (NEMO) is a critical constituent of the IKK complex (IKK-a, IKK-B, and IKK- 
y) involved in downstream of RIG-I and MDAS signaling, which is necessary for I-IFN pro- 
duction. The cleavage of NEMO caused by NSP5 impedes the synthesis of I-IFN [21]. In addi- 
tion to nonstructural proteins, accessory proteins such as NS6 encoded by PDCoV also play an 
important role in the inhibition of IFN-£ production by interfering with the binding of RIG-I/ 
MDAS5 and PDCoV [22]. However, the alteration of host cell innate immune response affected 
by PDCoV infection still needs to be explored. In particular, the main innate immune-associ- 
ated genes and signaling pathways of host cell modulated by PDCoV have not yet been uncov- 
ered, which can hinder further discoveries regarding the mechanism of PDCoV infection. 

MicroRNAs (miRNAs) are a group of short noncoding RNAs known as gene regulators 
that can modulate extensive bioprocesses by repressing the expression of their target genes at 
the translation level and/or mRNA level. It has been reported that the expression of miRNAs 
in host cells can be changed by multiple coronaviruses, such as severe acute respiratory syn- 
drome coronavirus (SARS-CoV), hepatitis C virus (HCV), porcine epidemic diarrhea virus 
(PEDV) and transmissible gastroenteritis virus (TGEV) [23, 24, 25, 26]. Some of the host cell 
endogenous miRNAs have been found to play important roles in the modulation of innate 
immune response affected by coronavirus. On one hand, host cell-derived miRNAs can 
enhance the innate immune response to inhibit the replication of coronaviruses. For excample, 
miRNA-221-5p (miR-221-5p) activates NF-«B signaling pathway which leads to the upregula- 
tion of I-IFN and the downstream ISGs as well as the suppression of PEDV replication [27]. 
On the other hand, cell endogenous miRNAs also can inhibit innate immune repsonse result- 
ing in coronaviruses innate immune evasion. For example, inositol-requiring enzyme 1 a 
(IRE1a:) modulates the decreasing of miR-30a-5p which involved in the signaling pathway of 
miR-30a-5p/SOCS1/3, resulting in the promotion of TGEV’s replication by escaping innate 
immune response [28]. Interferon stimulatory gene (IFITM1) is a target of miR-130a. Knock- 
down of miR-130a can increase the expression of IFITM1 resulting in the inhibition of HCV 
replication, which indicate that miR-130a directly involves in innate immune evasion of HCV 
[29]. As of yet, however, there has been no research into the host cell endogenous miRNAs 
which may involve in innate immune response modulation affected by PDCoV infection. 

To explore the modulation of innate immune response of host cell affected by PDCoV, this 
study investigated the global gene expression profiles of PK-15 cells at 0, 24, 36 h post infection 
(hpi) via RNA sequencing. Subsequently, differentially expressed genes screening, gene func- 
tion and signaling pathway enrichment analysis, and the prediction of PK-15 cell endogenous 
miRNAs that may make constributions to the innate immune response modulation affected 
by PDCoV were further performed. Our study provides a foundation for future systematic 
explorations of the PDCoV infection mechanism. 


Materials and methods 
PK-15 cell culture 


PK-15 cells, which were stored in our cell culture lab in this institute (Tianjin Institute of Ani- 
mal Husbandry and Veterinary Medicine, Tianjin Academy of Agricultural Sciences, Tianjin, 
China) and used in our previous research project [30], were used in this study. The PK-15 cells 
were cultured in Modified Eagle’s Medium (MEM) (HyClone) supplemented with 8% fetal 
bovine serum (ExCell), and 1% penicillin and streptomycin. Cell culture flasks were held at 
37°C ina humidified incubator set at 5% CO>. 
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PDCoV preparation and inoculation 


The PDCoV strain in this study was previously isolated by Swine Disease Research Innovation 
Lab of our Institute from small intestine specimens of piglets showed diarrheal syndromes 
[30]. Those piglets were raised in a swine farm in a county under Tianjin municipal adminis- 
tration. This PDCoV strain was propagated in PK-15 cell cultures. MEM supplemented with 
1% pancreatin (Sigma-Aldrich) was prepared for PDCoV dilution and cell culture after viral 
adsorption. PDCoV with a final concentration of 0.1 multiplicity of infection (MOI) was used 
for adsorption in single-line PK-15 cells for 2 h, followed by replacement with new cell culture 
medium. At 0, 24, and 36 hpi, cells were collected for RNA isolation. 


Indirect immunofluorescence assay (IFA) 


PDCoV infected PK-15 cells with 0.01 multiplicity of infection (MOI) and uninfected PK-15 
cells in 96-well plates were washed twice with phosphate-buffered saline (PBS), then fixed with 
100% ethanol for 30 min. Cells were incubated with rabbit antiserum (against PDCoV N pro- 
tein) at 1:200 dilution and the reaction was kept at a 37°C incubation for 1 h followed by 6 
times of washing with PBS. Then a secondary antibody of FITC-conjugated goat anti-rabbit 
antiserum (Solarbio, beijing, China) at 1:50 dilution was added and incubated at 37°C for 1 h. 
Cell wash was done in the same fashion as above. The IFA staining result was examined by 
using an fluorescent microscope. 


RNA isolation and quality control 


Total RNA was extracted using a GenElute™ Total RNA Purification Kit (Sigma-Aldrich) fol- 
lowing the manufacturer’s instructions. Polyacrylic acid amine gel electrophoresis was used to 
identify the RNA samples having no degradation or DNA pollution. Samples with OD (260/ 
280) ratios in the range of 1.8-2.0 and OD (260/230) ratios from 1.8 to 2.2 met the requirement 
of sequencing as identified using a NanoPhotometer spectrophotometer. RNA samples with 
RNA integrity numbers (RINs) greater than 7 were selected for subsequent RNA sequencing 
which was performed using an Agilent 2100 bioanalyzer. 


qRT-PCR 


Relative expression levels of 5 randomly selected genes were examined among the PK-15 cells 
at 0, 24, and 36 hpi by qRT-PCR using TB Green™ Premix Ex Taq™™ II (Tli RNaseH Plus) 
(Takara Biomedical Technology). The relative gene expression levels of the 5 genes were nor- 
malized to GAPDH expression using the 2-4CT method. The primers used in qRT-PCR were 
listed in S1 Table. The (RT-PCR reactions were performed on an ABI7500 StepOnePlus Real- 
Time PCR System (Thermo Fisher Scientific). 


cDNA library preparation and sequencing 


The mRNAs with poly(A) in 3’UTR were enriched from the total RNA using Magnetic Beads 
Oligo(dT) (Dynal Biotech ASA), then disrupted into short fragments using a Magnesium RNA 
Fragmentation Module kit (NEB). With fragmented mRNAs as templates and random oligo- 
nucleotides as primers, the first strands of cDNA were synthesized. Subsequently, the second 
strands of cDNA were generated using DNA polymerase I and deoxyribonucleotide triphos- 
phates (dNTPs). The terminals of purified double-strand cDNA sequences were repaired and 
provided with poly(A) and sequencing adapters. After ligating the adapters, the cDNA with 
200 bp were screened using AMPure XP beads. Polymerase chain reaction (PCR) was used for 
cDNA amplification and carried out on a Veriti™ 96-well Thermal Cycler (Thermo Fisher 
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Scientific). After further purifying the cDNA production using AMPure XP beads, a cDNA 
library was created. The insert sizes of cDNAs were detected using an Agilent 2100 bioanaly- 
zer. High-quality cDNAs with insert sizes ranging from 250 to 300 bp were selected for further 
quantification via (RT-PCR. According to the requirements of deep sequencing, cDNAs hav- 
ing target concentrations and data sizes were pooled, followed by RNA sequencing, which was 
performed using an Illumina HiSeqTM 4000 platform (Novogene Corp.). 


Clean reads filtering and quality control 


High-quality clean reads were screened from raw reads by removing reads with sequencing 
adapters or unknown nucleobases as well as low-quality reads in which the number of nucleo- 
bases with Qphrea quality scores < 20 were > 50% of the total read length. The distribution 
sequencing error rate was used to identify the quality of the data from Illumina sequencing. 
Q20 and Q30 represented the percent of nucleobases having 1% and 0.1% error rates, respec- 
tively, of the global nucleobase numbers determined using the Illumina sequencing platform. 


Genome mapping 


The swine reference genome (NCBI genome database: Sus scrofa 11.1) was used as the refer- 
ence genome for our RNA sequencing data. Genome mapping was performed using HISAT 
software (http://www.ccb.jhu.edu/software/hisat/index.shtml). The instructions of the HISAT 
software website was followed to prepare the data to be uploaded and also for the running 
analysis protocol. 


Differential expressed genes screening 


The DESeq2 software, available from the link at https://bioconductor.org/packages/release/ 
bioc/vignettes/DESeq2/inst/doc/DESeq2.html, was used for analysis of the differential expres- 
sion levels of genes in each comparison group based on Anders and Huber’s report [31]. The 
thresholds of the p-values were determined using the false discovery rate (FDR, also referred to 
aS Paqj). The threshold for differentially expressed gene screening was established to be a value 
of |log2(fold change in a comparison group)| > 1 and a false discovery rate (FDR, or pagj) < 
0.05. The FPKM (expected number of Fragments Per Kilobase of transcript sequence per Mil- 
lion base pairs sequenced) of each gene was calculated and used to represent the gene expres- 
sion value [32]. 


GO function and KEGG pathway enrichment analysis 


Function classification and enrichment analysis of DEGs were performed using GOseq soft- 
ware based on GO (http://geneontology.org/page/go-enrichment-analysis). Pathway enrich- 
ment was analyzed based on the KEGG database (https://www.kegg.jp/kegg/kegg1.html) and 
carried out using the KOBAS 3.0 web server. GO terms or KEGG pathways with p-values < 
0.05 were considered to be significantly enriched. 


Host cell endogenous miRNA prediction 


Host cell endogenous miRNAs were predicted using RNAhybrid (https://bibiserv.cebitec.uni- 
bielefeld.de/rnahybrid/), PITA (http://genie.weizmann.ac.il/pubs/mir07/mir07_prediction. 
html), and miRanda (http://www.microrna.org/microrna/home.do). The miRNAs that had 
reliable binding sites to target genes and commonly existed in 3 prediction results were consid- 
ered to be the contributors to the alteration of innate immune response induced by PDCoV. 
The binding sites between miRNA and target genes that had minimum free energy (MFE) 
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levels < —10 kcal/mol and p-values < 0.05 as determined by RNAhybrid were considered to be 
reliable binding structures. 


Statistical analysis 


The Student’s t test was used for statistical analyses in this study and carried out via SPSS Sta- 
tistics 24.0.0 (IBM, USA), as described in a previous study [33]. Data were presented as the 
mean + standard deviation (SD) of 2 independent experiments (n = 2). Data with p < 0.05 
were deemed to be statistically significant. 


Results 


Identification of cytopathic effect and PDCoV N protein in PDCoV- 
inoculated PK-15 cells 


PK-15 cells inoculated with PDCoV at 0.1 MOI were subsequently developed cytopathic effects 
(CPE) within 36 hpi. PDCoV induced CPEs on PK-15 cells was aggressively increased along 
with the increasing of PDCoV infection time (Fig 1). PDCoV-infected PK-15 cells at 24 hpi 
and 36 hpi yielded a same virus titer of 10*° TCIDso/mL. The IFA stained positive to PDCoV 
tested on these CPE cells (Fig 2). 


Overview of RNA sequencing results 


In order to further explore the alteration in the transcriptome profile of PK-15 cells during 
PDCoV infection, the global gene expression profiles of PK-15 cells after infection, including 
24 hpi (INF_24h) and 36 hpi (INF_36h), as well as uninfected PK-15 cells (NC), which acted 
as negative controls, were obtained by RNA sequencing. A total of approximately 1.6 billion 
raw reads were initially generated through sequencing. After filtering out the reads with adapt- 
ers, the reads with ambiguous nucleotides, and the reads of low quality, a total of nearly 156.2 
million clean reads remained, with approximately 7.8 G bases per sample. Clean data were 
indicated by a low sequencing error rate (0.03%) and the high quality of identified nucleobases 
via sequencing (Q20 > 97%, Q30 > 92%). The guanine-cytosine (GC)-content of clean reads 
ranged from 52.94% to 54.2%. These high-quality clean reads were subsequently mapped to a 
reference genome (NCBI genome database: Sus scrofa 11.1) which exhibited total mapping 
rates > 94% and unique mapping rates > 91%. From this, approximately 21,416 expressed 
genes, including 5,874 novel genes, were identified per sample. The RNA sequencing statistics 
are listed in Table 1. 


Fig 1. Cytopathic effect of PK-15 cell infections by PDCovV at 0, 24 and 36 hpi. In comparison with normal PK-15 cells (A), PDCoV infected cells at 24 hpi (B) 
exhibited CPE of rounded cells (arrow), shrinkaged cells and released cells from monolyer; after prolonged to 36 hpi, PDCoV infected showed more progressive 
development of CPEs (C). 


https://doi.org/10.1371/journal.pone.0223177.g001 
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Fig 2. IFA staining on PK-15 cells. IFA stained positive (green fluorescence) on PDCoV-infected PK-15 cells (A) and 
negative on normal PK-15 cells (B). 


https://doi.org/10.1371/journal.pone.0223177.g002 


Screening of differentially expressed genes 


The expression of all identified genes obtained through RNA sequencing were respectively 
compared for the group of uninfected PK-15 cells and cells at 24 hpi (INF_24h versus NC) as 
well as the group of PK-15 cells at 24 and 36 hpi (INF_36h versus INF_24h). Genes which dif- 
ferentially expressed in each comparison group with |log2(fold change in a comparison 
group)| > 1 and false discovery rate (FDR, or Pag;) < 0.05 were considered to be differentially 
expressed genes (DEGs). For INF_24h versus NC, there were 3,762 DEGs, with 2,660 genes 
upregulated and 1,102 genes downregulated in INF_24h (Fig 3A; S2 Table). For INF_36h ver- 
sus INF_24h, 560 DEGs were screened. Compared to INF_24h, 156 and 404 genes were 
respectively upregulated and downregulated in INF_36h (Fig 3B; $3 Table). The DEGs in these 
2 comparisons were the genes that were mainly affected during PDCoV infection. 

In order to identify the alteration of innate immune-related genes in host cells as PDCoV 
post infection time increased, porcine innate immune genes were further screened. Up to date, 
porcine innate immune genes have not been collected into a database which can be conve- 
niently updated frequently. In addition, the anatomy, organization, and response of pig 
immune system are very similar to humans’ [34]. Thus, the porcine innate immune genes that 
included in the DEGs were screened by referring to the innate immune genes list of the Inna- 
teDB (http://www.innatedb.ca/redirect.do?go=aboutIDB) which provides a total of 1,040 
curated human genes related to the innate immune response. There were 156 and 23 swine 
innate immune genes that were differentially expressed in INF_24h versus NC and INF_36h 
versus INF_24h, respectively. In INF_24h versus NC, innate immune DEGs with 120 upregu- 
lated and 36 downregulated tendencies were identified in INF_24h compared to NC (Fig 4A). 
In INF_36h versus INF_24h, 19 genes were downregulated in INF_36h while 4 genes were 
upregulated (Fig 4B). Of note, most of the differential expressed innate immune genes includ- 
ing I-IFNs such as IFN-a@1 and IFN-1 (also known as IFN-f) were upregulated in cells at 24 
hpi comparing to uninfected cells, and downregulated in cells at 36 hpi comparing to cells at 


Table 1. Statistics of RNA sequencing results. 


Sample | Raw reads Clean Clean Sequencing error | Q20 | Q30 | GCcontent | Total mapping | Unique mapping | Expressed Novel 
reads bases (G) rate (%) (%) (%) (%) rate (%) rate (%) genes genes 
NC_1 26,873,398 | 26,397,492 7.92 0.03 97.18 | 92.50 53.06 95.18 92.56 21516 5591 
NC_2 24,489,789 | 23,998,778 7.20 0.03 97.19 | 92.54 52.94 94.80 92.25 20813 5386 
INF_24h1 | 29,972,480 | 29,121,085 8.74 0.03 97.46 | 93.44 53.61 94.10 91.20 22169 6160 
INF_24h2 | 25,687,216 | 24,943,590 7.48 0.03 97.58 | 93.73 53.85 94.39 91.53 21425 5998 
INF_36h1 | 29,874,857 | 28,998,514 8.70 0.03 97.60 | 93.72 54.20 94.32 91.54 21349 6097 
INF_36h2 | 23,435,580 | 22,780,078 6.83 0.03 97.16 | 92.53 53.56 94.50 91.71 21224 6012 


https://doi.org/10.1371/journal.pone.0223177 001 
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Fig 3. DEGs in INF_24 h versus NC and INF_36 h versus INF_24 h. (A): DEGs upregulated and downregulated on 
PK-15 cells at 24 hpi represented by red dots, uninfected PK-15 cells represented by green dots. Genes with no 
significant differences in expression represented in blue. (B): DEGs upregulated in INF_36 h compared to INF_24h, 
the same color coded dots as (A). 


https://doi.org/10.1371/journal.pone.0223177.g003 


24 hpi, which indicated that innate immune response was initially enhanced and then inhib- 
ited in cells at 24 hpi and 36 hpi, respectively. 

In addition, these 2 comparison groups shared 400 DEGs (Fig 5A). A total of 17 DEGs includ- 
ing IFN-a1 and IFN-1 belonged to porcine innate immune-associated genes (Fig 5B). Five of 
the 17 genes were randomly selected for validation via qRT-PCR, and the results showed that the 
relative expressions of these genes were coincident with our RNA sequencing results (Fig 6). 


GO function enrichment analysis 


The DEGs in the 2 comparison groups were annotated by the Gene Ontology (GO) function 
database, which divided them into 3 macroscopic groups, including biological process (BP), 
cellular component (CC), and molecular function (MF). The number of DEGs in the INF_24h 
versus NC comparison that fell into each Gene function classification were higher than the 
number of DEGs in the INF_36h versus INF_24h comparison (Fig 7). Also, in relation to 
INF_36h versus INF_24h, the DEGs annotated by GO terms in INF_24h versus NC were 
more abundant (Fig 8). 

To explore the function of genes that affected during the alteration of innate immune 
response caused by PDCoV, DEGs whose relative expression tendencies were consistent with 
most of the differential expressed host innate immune genes screened in the present study 
were further analyzed using GO function enrichment. Namely, DEGs upregulated in cells at 
24 hpi compared to uninfected cells and DEGs downregulated in cells at 36 hpi compared to 
cells at 24 hpi were selected for further analysis with GO function enrichment. For INF_24h 
versus NC, genes which were upregulated in cells at 24 hpi were enriched in 99 terms of BP, 3 
terms of CC, and 9 terms of MF (S4 Table). Most of these genes were associated with response 
against viruses, host defense response, and immune response (Fig 8A). In another comparison 
group, genes with downregulated tendencies in cells at 36 hpi were enriched in 65 terms of BP, 
1 term of CC, and 2 terms of MF (S5 Table). These genes mainly enriched in the functions 
related to host defense responses, such as defense response, response to external biotic stimuli, 
response to external stimuli, and in the functions of cytokines involved in immune and innate 
immune response, such as cytokine activity and cytokine receptor binding (Fig 8B). Since 
genes annotated with GO terms associated with immune and innate immune responses were 
commonly shared, the enriched immune and innate immune response-related GO terms were 
both selected to represented (Tables 2 and 3). 
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Fig 4. Innate immune-associated DEG. Innate immune-associated genes of DEGs in INF_24h versus NC (A) and 
INF_36h versus INF_24h (B) were illustrated in heatmaps that were created based on the expression of each DEG. (A) 
Genes with upregulated tendencies in INF_24h comparing to NC were displayed in red while genes with 
downregulated tendencies in INF_24h comparing to NC were shown in blue. (B) Genes with upregulated tendencies 


in INF_36h comparing to INF_24h were represented in red while genes with downregulated tendencies were shown in 
blue. 


https://doi.org/10.1371/journal.pone.0223177.g004 


KEGG pathway enrichment analysis 


DEGs used in GO function enrichment analysis were also selected for analysis in Kyoto Ency- 
clopedia of Genes and Genomes (KEGG) pathway enrichment. A total of 36 signaling path- 
ways (Fig 9A), including 6 innate immune response-associated pathways (Table 4), were 
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Fig 5. Common genes belonging to DEGs in INF_24h versus NC and INF_36h versus INF_24h. (A) A total of 400 
DEGs were commonly expressed in the 2 comparison groups. (B) Innate immune-associated genes of DEGs commonly 
shared in cells at 0, 24, and 36 hpi were illustrated in a heatmap. The heatmap was constructed based on the expression 
for each DEG with 2 replicates examined. 


https://doi.org/10.1371/journal.pone.0223177.g005 


enriched in INF_24h versus NC. Additionally, a total of 26 signaling pathways (Fig 9B), includ- 
ing 9 related to innate immune response (Table 5), were enriched in INF_36h versus INF_24h. 
This finding indicates that extensive bioprocesses were affected by PDCoV infection. Of note, 
5 signaling pathways were commonly represented in the 2 enrichment results. These innate 
immune-associated pathways were mainly affected during PDCoV enhanced or inhibited 
innate immune response of PK-15 cells. In addition, 3 of the 5 pathways belong to pattern rec- 
ognition receptors (PRRs) mediated signalings, including NOD-like receptor signaling path- 
way, RIG-I-like receptor signaling pathway, and cytosolic DNA-sensing pathway, which can 
mediate the J-IFN production. One of the 5 pathways, Jak/STAT signaling pathway, belongs to 
the downstream signaling of I-IFN response which can mediate the production of ISGs. These 
4 pathways are known to play important roles in the process that I-IFN develops antiviral func- 
tion in innate immune response. 
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Fig 6. Relative expression of 5 randomly selected genes among PK-15 cells at 0, 24 and 36 hpi. Relative expression of all of the genes were 
coincident with our RNA sequencing results. Data are displayed as mean + SD (n = 2); *p < 0.05; **p < 0.01; ***p < 0.001. 


https://doi.org/10.1371/journal.pone.0223177.g006 


Host cell endogenous miRNA prediction 


It is well known that host cell endogenous miRNAs are crucial contributors in modulating 
host cells innate immune response which affected by coronaviruses. Given this, the present 
study further predicted PK-15 cell-derived miRNAs which may involved in the alteration of 
PK-15 cell innate immune response caused by PDCoV infection. Considering the importance 
of the 4 innate immune signaling pathways in the process that I-IFN developed antiviral func- 
tion, DEGs which were covered in these signaling pathways were collected for miRNA predic- 
tion. A total of 16 miRNAs were likely corelated with DEGs, of which 12 could modulate 1 or 
more innate immune signaling pathway by directly regulating unique target genes (Table 6), 
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Fig 7. Gene function classifications of DEGs in 2 comparison groups. Gene function classifications of DEGs in INF_24h versus NC and INF_36h versus 


INF_24h comparisons are represented by blue and red bars, respectively. 


https://doi.org/10.1371/journal.pone.0223177.g007 
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Fig 8. GO function enrichment results of DEGs in 2 comparison groups. Top 10 significantly enriched GO terms of each category, including biological 
process (BP), cellular component (CC), and molecular function (MF) were displayed in this figure (A: INF_24h versus NC; B: INF_36h versus INF_24h). 
Enriched GO terms with p < 0.05 were considered to be statistically significant. 


https://doi.org/10.1371/journal.pone.0223177.g008 


and 4 could modulate 2 or more innate immune signaling pathways by directly regulating 
multiple target genes (Table 7; Fig 10). 


Discussion 


TGEV, PEDV, and PDCoV are coronaviruses causing most severe swine enteric disease prob- 
lems and resulting in significant economic losses in swine industry [12]. It has been reported 
that pigs experiencing severe diarrhea are frequently coinfected with TGEV, PEDV, and 
PDCoV. However, even in the absence of TGEV and PEDV, PDCoV has high pathogenicity, 
making it worthy of significant concern. In general, PDCoV was the only pathogenic factor 
among these three viruses, which could be detected in feces and intestinal samples and associ- 
ated with outbreaks of swine diarrheal disease in the U.S. in 2014 [12]. To date, the infection 
mechanisms of TGEV and PEDV at the mRNA or miRNA levels have already been revealed 
by high-throughput sequencing, which has delineated the key genes (or miRNAs) and signal- 
ing pathways involved in the process of viral infection [26, 35]. However, analogous PDCoV 
research still lags behind. The present study is the first related to the discovery of the infection 
mechanism of PDCoV at the transcriptome level by RNA sequencing. 

Escaping the host cell innate immune response is the key stage in the successful infection of 
PDCoV in host cells. To date, research studies have made great progress revealing the strate- 
gies that PDCoV utilizes to escape I-IFN-associated host cell innate response. Nevertheless, we 
still need further research studies to completely understand the complex molecular mechanism 
of the alteration of innate immune response affected by PDCoV. Given this, our study focused 
primarily on discovering innate immune-associated genes and signaling pathways that altered 
by PDCoV infection in PK-15 cells. 
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Table 2. GO function enrichment results of upregulated DEGs in cells at 24 hpi comparing to unfected cells. 


Category ID Description P-value Gene name Count 
BP GO:0002252 Immune effector process 1.47E-05 | RSAD2, OAS1, DDX58, MX2, STAT1, STAT2, C1S, IL27, IFNA1, IRF1, TRAF6, 15 
ISG20, IL10, AZU1 
BP GO:0002376 Immune system process 2.42E-05 | RSAD2, OAS1, DDX58, MX2, STAT1, STAT2, OLR1, PSMB8, SOX9, CCL5, SDC4, 31 
C1S, CCL2, IL15, IL27, ILI2A, IFNA1, IRF1, TRAF6, ISG20, IL10, AZU1, ESR1 
BP GO:0006955 Immune response 8.04E-05 | RSAD2, OAS1, DDX58, MX2, STAT1, STAT2, CCL5, C1S, CCL2, IL15, IL27, IL12A, 25 
IRF1, TRAF6, ISG20, IL10, AZU1, ESR1 
BP GO:0001819 | Positive regulation of cytokine 1.45E-03 | RSAD2, DDX58, STAT1, IL27, IRF1, TRAF6, IL10, AZU1 8 
production 
BP GO:0045087 Innate immune response 1.70E-03 | RSAD2, OAS1, DDX58, MX2, STAT1, STAT2, CCL5, C1S, CCL2, IL27, IRF1, ISG20, 13 
ESRI 
BP GO:0032635 Interleukin-6 production 2.99E-03 | DDX58, GHSR, TRAF6, NOS2, IL10 5 
BP GO:0001817 Regulation of cytokine 4.50E-03 | RSAD2, DDX58, STAT 1, IL27, GHSR, IRF1, TRAF6, NOS2, IL10, AZU1 10 
production 
BP GO:0042035 Regulation of cytokine 4.97E-03 | IL27, GHSR, IRF1, TRAF6, AZU1 5 
biosynthetic process 
BP GO:0042089 Cytokine biosynthetic process 7.74E-03 | IL27, GHSR, IRF1, TRAF6, AZU1 5 
BP GO:0042107 Cytokine metabolic process 7.74E-03 | IL27, GHSR, IRF1, TRAF6, AZU1 5 
BP GO:0001816 Cytokine production 1.24E-02 | RSAD2, DDX58, STAT1, IL27, GHSR, IRF1, TRAF6, NOS2, IL10, AZU1 10 
BP GO:0002758 Innate immune response- 2.05E-02 | RSAD2, DDX58, IRF1, ESR1 4 
activating signal transduction 
BP GO:0002218 Activation of innate immune 2.93E-02 | RSAD2, DDX58, IRF1, ESR1 4 
response 
BP GO:0007249 I-kappaB kinase/NF-kappaB 2.93E-02 | STAT1, OPTN, TRAF6, ESR1 4 
signaling 
BP GO:0045089 Positive regulation of innate 2.93E-02 | RSAD2, DDX58, IRF1, ESR1 4 
immune response 
BP GO:0002682 Regulation of immune system 2.99E-02 | RSAD2, DDX58, STAT 1, SOX9, CCL5, SDC4, C18, CCL2, IL27, IRF1, TRAF6, IL10, 13 
process ESRI 
BP GO:0002757 Immune response-activating 3.74E-02 | RSAD2, DDX58, IRF1, TRAF6, ESR1 5 
signal transduction 
BP GO:0002697 | Regulation of immune effector | 4.71E-02 | RSAD2, DDX58, IL27, TRAF6, IL10 5 
process 
MF GO:0005126 Cytokine receptor binding 9.74E-04 | CCL5, CCL2, IL15, IL27, BDNF, IL12A, IFNA1, KITLG, TRAF6 13 
MF GO:0005125 Cytokine activity 1.86E-02 | NAMPT, CCL5, CCL2, IL15, IL27, IL12A, IFNA1, IL10 10 


Immune and innate immune response-related GO terms were represented in this table. 


https://doi.org/10.1371/journal.pone.02231 77.1002 


Reported studies indicated that multiple cell lines including ST, LLC-PK and IPI-2I were 
permissive to PDCoV [36,37]. We found the PK-15 cell line was permissive to PDCoV in this 
study. However, although the degree of CPE in PK-15 cells from 24 hpi to 36 hpi showed an 
increasing tendency, cells at 12 h different post infection (pi) time yielded a similar virus titer. 
The irrelevance between the number of infected cells and virus titer at each pi time point could 
result from susceptibility of PK-15 to PDCoV, which was different from other cell lines of per- 
missive to PDCoV. 

To further investigate the mechanism of PDCoV infection, all of the expressed genes were 
identified in PK-15 cells at 0, 24 and 36 hpi through RNA sequencing. By comparing the 
expressed genes in cells at 0 and 24 hpi, we found that most of the innate immune-associated 
genes were significantly upregulated at 24 hpi comparing to cells without PDCoV infection. 
The results indicated that I-IFN-mediated innate immune response was triggered during 
PDCoV at an early infection time point (24 hpi). In the comparison between cells at 24 hpi 
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Table 3. GO function enrichment results of downregulated DEGs in cells at 36 hpi comparing to cells at 24 hpi. 


Category ID Description P-value Gene name Count 
BP GO:0002376 Immune system process 1.69E-02 CCL2, CCL5, IFNA1, AZUI1, IL27 7 
BP GO:0006955 Immune response 1.87E-02 CCL2, CCL5, AZU1, IL27 6 
BP GO:0042035 Regulation of cytokine biosynthetic process 2.19E-02 AZUI, IL27 2 
BP GO:0042089 Cytokine biosynthetic process 2.59E-02 AZUI, IL27 2 
BP GO:0042107 Cytokine metabolic process 2.59E-02 AZUI, IL27 2 
BP GO:0001817 Regulation of cytokine production 3.74E-02 AZU1, IL27, NOS2 3 
MF GO:0005125 Cytokine activity 2.78E-03 CCL2, CCL5, IFNA1, IL27 5 
MF GO:0005126 Cytokine receptor binding 3.50E-03 CCL2, CCL5, IFNA1, IL27 5 


Immune and innate immune response- 


related GO terms were represented in this table. 


https://doi.org/10.1371/journal.pone.02231 77.1003 


and 36 hpi, the PRRs and ISGs described above which were elevated significantly in INF_24h 
versus NC were not downregulated significantly in cells at 36 hpi comparing to cells at 24 hpi, 
but the significant inhibition of I-IFNs expression still can be found. The result indicated that 
the I-IFNs production was reduced at 36 hpi comparing to 24 hpi which also indicated that the 
I-IFN response might be inhibited at 36 hpi induced by PDCoV. The reason that the PRRs or 
ISGs described above were not significantly down-regulated at this infection time point may 
due to there are other key molecules which regulated by PDCoV for suppression of I-IFNs pro- 
duction or I-IFN response at 36 hpi. Besides, the results also indicated that the genes which 
regulated by PDCoV in a relative early infection time (24 hpi) and a relative late infection time 


(36 hpi) for innate immune response modulation were different. 


The results of GO function and KEGG signaling pathway enrichment of differentially 
expressed genes possessing the expression tendency consistent with most of the innate 
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Fig 9. KEGG pathway enrichment results of DEGs in 2 comparison groups. KEGG signaling pathways which were significantly enriched in 
INF_24h versus NC and INF_36h versus INF_24h are represented in (A) and (B), respectively. Enriched signaling pathways with p < 0.05 were 
considered to be statistically significant. 
https://doi.org/10.1371/journal.pone.0223177.g009 
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Table 4. KEGG pathway enrichment results of upregulated DEGs in INF_24h comparing to NC. 


ID Description Gene | P-value Gene name Count 
ratio 
ssc04060 Cytokine-cytokine 27/434 | 3.46E- | IFNB1, CXCL10, IFN-ALPHAOMEGA, TNFRSF12A, CXCL11, CCL5, TNFSF13B, TNFRSF1B, 27 
receptor interaction 04 CCL2, IL15, CCR7, FLT4, IL22RA1, IL12A, IFNA1, KITLG, IL10RA, IL10, IL12RB1, AMHR2, 
CXCL12, TNFRSF9, IL11, HGF, MPL, AMH, TNFRSF11A 
ssc04621 NOD-like receptor 21/434 | 1.01E- | GBP1, OASI, IFNB1, STAT1, STAT2, IFN-ALPHAOMEGA, LOC100155195, TXNIP, CASP1, 21 
signaling pathway 03 NAMPT, LOC100156073, IRF7, CCL5, LOC100523492, CCL2, TNFAIP3, IFNA1, TRAF6, BIRC3, 
ERBIN, XIAP 
ssc04064 | NF-kappa B signaling 13/434 | 1.76E- | DDX58, TRIM25, PTGS2, TNFSF13B, VCAM1, CARD11, TNFAIP3, CFLAR, TRAF6, BIRC3, 13 
pathway 03 CXCL12, XIAP, TNFRSF11A 
ssc04622 RIG-I-like receptor 11/434 | 4.75E- | IFNB1, DDX58, CXCL10, IFN-ALPHAOMEGA, TRIM25, CASP10, IRF7, CYLD, IL12A, IFNA1, 11 
signaling pathway 03 TRAF6 
ssc04630 Jak-STAT signaling 19/434 | 6.74E- | IFNB1, STAT1, STAT2, IFN-ALPHAOMEGA, SOCS1, FHL1, IL29, PIM1, IL15, IL22RA1, IL12A, 19 
pathway 03 IFNAI, CISH, ILI0RA, IL10, ILI2RB1, STAT4, IL11, MPL 
ssc04623 | Cytosolic DNA-sensing 8/434 | 2.43E- | IFNB1, DDX58, CXCL10, IFN-ALPHAOMEGA, CASP1, IRF7, CCL5, IFNA1 8 
pathway 02 


Innate immune response-related pathways were represented in this table. 


https://doi.org/10.1371/journal.pone.02231 77.1004 


immune-associated genes in INF_24h versus NC as well as INF_36h versus INF_24h revealed 
that extensive gene functions and signaling pathways including innate immune-related ones 
were affected at different infection time of PDCoV. Additionally, we found that the GO terms 
and KEGG signaling pathways which enriched in two comparison groups were distinct. These 
findings suggested that the gene functions and signaling pathways which affected by PDCoV 
for modulation of innate immune response at different post infection times were complicated 
and different. 
Five innate immune-associated signaling pathways which were commonly enriched in 2 

comparison groups are primarily affected during the alteration of innate immune response of 


Table 5. KEGG pathway enrichment results of downregulated DEGs in INF_36h comparing to INF_24h. 


ID Description Gene P-value Gene name Count 
ratio 

ssc04620 Toll-like receptor signaling 8/75 1.77E- | CXCL10, IFN-ALPHAOMEGA, IFNB1, CCL5, CXCL11, IFNA1, RAC1, MAP2K6 8 
pathway 05 

ssc04623 | Cytosolic DNA-sensing pathway 5/75 3.28E- | CXCL10, IFN-ALPHAOMEGA, IFNB1, CCL5, IFNA1 5 
04 

ssc04060 Cytokine-cytokine receptor 8/75 2.53E- | CXCL10, IFN-ALPHAOMEGA, IFNB1, CCL2, CCL5, CXCL11, IFNA1, CCR7 8 
interaction 03 

ssc04151 PI3K-Akt signaling pathway 10/75 4.74E- | IFN-ALPHAOMEGA, IFNB1, YWHAZ, IFNA1, RAC1, FGF18, GNG11, PDGFC, 10 

03 COL9A2, COL4A6 

ssc04622 RIG-I-like receptor signaling 4/75 7.52E- | CXCL10, IEFN-ALPHAOMEGA, IFNB1, IFNA1 4 
pathway 03 

ssc04015 Rapl signaling pathway 7175 9.37E- | RALGDS, CTNND1, SIPAIL2, RACI1, FGF18, MAP2K6, PDGFC 7 
03 

ssc04630 Jak-STAT signaling pathway 5/75 3.14E- | IFN-ALPHAOMEGA, IFNB1, IL29, FHL1, IFNA1 5 
02 

ssc04621 NOD-like receptor signaling 5/75 3.35E- | IFN-ALPHAOMEGA, IFNB1, CCL2, CCL5, IFNA1 5 
pathway 02 

ssc04014 Ras signaling pathway 6/75 4.82E- | RALGDS, RAC1, REL, FGF18, GNG11, PDGFC 6 
02 


Innate immune response-related pathways were represented in this table. 


https://doi.org/10.1371/journal.pone.0223177.t005 
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Table 6. Predicted miRNAs with unique targets. 


Species MiRNA Target gene Innate immune signaling pathway 

Sus scrofa miR-195 IFN-ALPHAOMEGA NOD-like receptor signaling pathway, RIG-I-like receptor signaling pathway, 
Jak-STAT signaling pathway, Cytosolic DNA-sensing pathway 

Sus scrofa miR-7137-3p OAS1 NOD-like receptor signaling pathway 

Sus scrofa miR-4334-5p TNFAIP3 NOD-like receptor signaling pathway 

Sus scrofa miR-139-3p TNFAIP3 NOD-like receptor signaling pathway 

Sus scrofa miR-493-3p IFNA1 NOD-like receptor signaling pathway, RIG-I-like receptor signaling pathway, 
Jak-STAT signaling pathway, Cytosolic DNA-sensing pathway 

Sus scrofa miR-9822-3p CASP1 NOD-like receptor signaling pathway, Cytosolic DNA-sensing pathway 

Sus scrofa let-7a CISH Jak-STAT signaling pathway 

Sus scrofa miR-4335 TNFAIP3 NOD-like receptor signaling pathway 

Sus scrofa miR-1306-5p TRIM25 RIG-I-like receptor signaling pathway 

Sus scrofa miR-9853-5p IL12RB1 Jak-STAT signaling pathway 

Sus scrofa miR-92b-5p TNFAIP3 NOD-like receptor signaling pathway 

Sus scrofa miR-671-5p CYLD RIG-I-like receptor signaling pathway 


https://doi.org/10.1371/journal.pone.02231 77.1006 


host cells which caused by PDCoV infection. Of note, 4 of the 5 signaling pathways including 
NOD-like receptor signaling pathway, RIG-I-like receptor signaling pathway, Jak-STAT sig- 
naling pathway, and cytosolic DNA-sensing pathway play crucial roles in the process that 
I-IFN develops antiviral function in innate immune response. Moreover, NF-kappa B signaling 
pathway is specifically enriched in INF_24h versus NC, while Toll-like receptor signaling path- 
way, PI3K-Akt signaling pathway, Rap] signaling pathway and Ras signaling pathway are spe- 
cifically enriched in INF_36h versus INF_24h. Recent studies have reported that PDCoV 
works against the RIG-I-like receptor signaling pathway and Jak/STAT signaling pathway, 
which are located upstream and downstream of I-IFN respectively, in order to escape the host 
cell innate immune response [20, 21, 22]. However, the other innate immune-related signaling 
pathways that are affected by PDCoV infection against the host cell innate immune response 
discovered in our present study also deserve attention and should be explored further in the 
future. 

Since host cell endogenous miRNAs are important contributors to the modulation of cell 
innate immunity by coronaviruses [23-29], porcine miRNAs belonging to PK-15 cells that 


Table 7. Predicted miRNAs with multiple targets. 


Species | MiRNA | Target gene Innate immune signaling pathway 


Sus scrofa | miR-17-3p CCL2 NOD-like receptor signaling pathway 

STAT2 | NOD-like receptor signaling pathway, Jak-STAT signaling pathway 
Sus scrofa | miR-328 OAS1 NOD-like receptor signaling pathway 

TRIM25 | RIG-I-like receptor signaling pathway 

Sus scrofa | miR-874 ILIORA _| Jak-STAT signaling pathway 


IFNA1 NOD.-like receptor signaling pathway, RIG-I-like receptor signaling pathway, 
Jak-STAT signaling pathway, Cytosolic DNA-sensing pathway 


ILI1 Jak-STAT signaling pathway 
Sus scrofa | miR-370 IL12RB1 | Jak-STAT signaling pathway 
FHL1 Jak-STAT signaling pathway 


IRF7 NOD-like receptor signaling pathway, RIG-I-like receptor signaling pathway, 
Cytosolic DNA-sensing pathway 


TXNIP | NOD-like receptor signaling pathway 


https://doi.org/10.1371/journal.pone.0223177 t007 
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position: 62 
target: CCL2 5' -UCUGCAGCCCUCCUGUGCCUGCUGCUCACUGCAGC- 3' MFE: -32.3 kcal/mol 
LH HILL LETT p value: 1.50E-02 
miRNA: miR-17-3p 3' -GAUGUU---—-------- CACGGA—-——— AGUGACGUCA- 5' 
position: 730 
target: STAT2 5' -UCAAGUGUUGUCACUGCAGC- 3° MFE: -30.8 kcal/mol 
HTL TV TT TEP LT EIA p value: 4.85E-02 
miRNA: miR-17-3p 3' -GAUGUUCACGG AAGUGACGUCA- 5' 
position: 715 
target: OAS1 5' -AGGGUCGGGGCCCGGAGCCUUCAUCCCUAGAGGGCCAU- 3' MFE: -33.1 kcal/mol 
HE UIE TL HLL p value: 3.56E-02 
miRNA: miR-328 3'-UGCCU----UCCCG----UCU-----------------------— CUCCCGGUC- 5' 
position: 835 
target:TRIM25 5' -UGCGGUGGGGGGGGGGGGGUAGGA- 3° MFE: -36.9 kcal/mol 
HL WEE LILLE p value: 6.70603 
miRNA: miR-328 3' -UGCCU-UCCC-GUCUCUCCCGGUC- 5° 
position: 201 


target: ILLORA 5' -GUGGUCAGUUGUGUGACCUUGGGCCAGCCUCUUCUCAGGGGGCC- 3' MFE:-34.5 kcal/mol 
l| | | | | | l| lI | | p value: 3.41E-02 


miRNA: miR-874 3' -CAGCCAG----—---------—----- GGAGCCCGGU--------------—----- CCCGUC- S' 

position: 521 

target: IFNA1 5' -CCUGUGCCUGGGAGAUCGUCAGGGCAGA- 3' MFE: -33.8 kcal/mol 

LETT p value: 2.41E-02 

miRNA: miR-874 3' -CAGCCAGGGAGCC---------—- CGGUCCCGUC- 5' 

position: 316 

target: IL-11 s' -~-GGGUCCCUGGGGGGCCUUGGGGGAGCAGGGCAGC- 3° MFE: -40.4 kcal/mol 

WET TH LH] p value: 2.88E-03 

miRNA: miR-874 3' -CAGCCAGGGAG----- CCOGG--—————-———— ucccGuc- 5° 

position: 726 

target: FHL1 5‘ -AGCGGGGAGGCCACGGUCCCGGCAGGCA- 3' MFE: -38.7 kcal/mol 

HETIL TEE p value: 2.87E-03 

miRNA: miR-370 3' -UGGUCCAA--GGUG----GGGUCGUCCG - _ 5' 

position: 1 

target: IL12RB1 5' -CAGGUUUGGCAGUGGACCAGCCCAGCAG-A- 3' MFE: -33.2 kcal/mol 
|| || I ITT p value: 3.29E-02 

miRNA: miR-370 3'-UGGUCCAA----------------- GGUGGGGUCGUCCG- 5' 

position: 416 

target: IRF7 5' -UCCGGGCCCGGCU GGCGCGGCUGCAGAGCCCCAGCAGGUA- 3' MFE: -36.3 kcal/mol 
HII lI HET TT p value: 8.24£-03 

miRNA: miR-370 3° -UGGUCCAA----- ‘GG--—- UGGGGUCGUCCG- 5° 

position: 676 

target: TXNIP 5' -GUGGGGUACAGAUCCACCCCAGC--C- 3 MFE: -33.1 kcal/mol 
| | || WT p value: 3.44E-02 

miRNA: ssc-miR-370 3" -UGGUCCA---------- AGGUGGGGUCGUCCG- 5S' 


Fig 10. Binding sites between miRNAs and their target genes covered in innate immune signaling pathways. Four 
miRNAs with multiple targets were selected to display the binding sites (A: miR-17-3p; B: miR-328; C: miR-874; D: 
miR-370). The results of binding sites were predicted by RNAhybrid. The binding sites with a minimum free energy 
(MFE) < —10 kcal/mol and a p-value < 0.05 were considered to be stable binding site structures. 


https://doi.org/10.1371/journal.pone.0223177.g010 


may modulate the 4 important PRRs-mediated innate immune signalings by targeting unique 
or multiple genes involved in these signalings were predicted in present study, which provides 
the foundation for a new research field concerning the investigation of the PDCoV infection 
mechanism. However, whether these miRNAs are required for PDCoV-induced innate 
immune modulation by regulating the potential target genes listed in the results (Tables 6 and 
7) still need to be further explore. 


Conclusions 


The preliminary results described in this report explore the innate immune-associated genes 
and signaling pathways affected by PDCoV infection, which provide the potential host cell 
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endogenous miRNAs that may contribute to the modulation of innate immune response 
induced by PDCoV infection through the regulation of innate immune genes and pathways. 
This study will lead to further research studies into the mechanism of PDCoV infection as well 
as the prevention of PDCoV-induced swine diarrheal disease. 
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